#Note: this code estimates Models 1-3 as reported in Table A.22 of the Supplemntal Appendix.
#   To do so, it reads-in the R data file "tvdata.rda" and then prints the relevant table output to screen.
#   Because these models take 1-2 weeks to run, the replication folder also includes a saved R workspace with
#   the final model output for the script below ("Table_A22_Models_1_to_3").
#   This saved R workspace can be loaded into R via load("Table_A22_Models_1_to_3") and once loaded the output
#   from the models below can be printed to screen via the summary commands that appear at the end of this script.

################################################################
#############################Set-up#############################
################################################################

#set memory size
memory.size(max=TRUE)
memory.limit(size=3000000)

#clear memory
rm(list=ls())

#load relevant packages                                             
library(foreign)
require(devtools)
install_version("amen", version = "1.3", repos = "http://cran.us.r-project.org")
library(amen)

#set seed
seed <- 12345
set.seed(seed)

#set working directory#
setwd("")

###########################################################
#################Load Dataset for Analysis#################
###########################################################

#load data
load('tvdata.rda')

############################################################
######################Conduct Analysis######################
############################################################

#set mcmc params
imps = 100000
brn = 50000
ods = 15
latDims = 2

#estimate Table A.22, Model 1
ameFit.1 = ame_rep(Y=Y, Xdyad=Xd,  Xrow=Rownode, Xcol=Colnode,  model='bin', symmetric=FALSE, R=latDims,
        nscan=imps, seed=seed, burn=brn, odens=ods, 
        rvar=TRUE,cvar=TRUE,dcor=FALSE, nvar=FALSE, plot=FALSE, print=FALSE) 
summary(ameFit.1)

#estimate Table A.22, Model 2
ameFit.2 = ame_rep(Y=Y, Xdyad=Xd,  Xrow=Rownode, Xcol=Colnode,  model='bin', symmetric=FALSE, R=latDims, 
        nscan=imps, seed=seed, burn=brn, odens=ods, 
        rvar=TRUE, cvar=TRUE, dcor=TRUE,nvar=FALSE, plot=FALSE, print=FALSE) 
summary(ameFit.2)

#estimate Table A.22, Model 3
ameFit.3 = ame_rep(Y=Y, Xdyad=Xd,  Xrow=Rownode, Xcol=Colnode,  model='bin', symmetric=FALSE, R=latDims, 
        nscan=imps, seed=seed, burn=brn, odens=ods, 
        rvar=TRUE, cvar=TRUE, dcor=TRUE,nvar=TRUE, plot=FALSE, print=FALSE) 
summary(ameFit.3)

#now summarize Models 1, 2 and 3 for Table A.22
summary(ameFit.1)
summary(ameFit.2)
summary(ameFit.3)
